module ModuleSolve_GLO_LOC
    
    use Globals
    use Toolbox
    use Functions_wa

contains

    subroutine Solve_GLO_LOC
    
    implicit none
    
    ! Local variable declarations
    real(8), dimension(NS) :: xra_hmax      ! maximum asset choice
    real(8) :: errVold                      ! old error in value function
    real(8) :: derrV                        ! change in error
    integer :: itV, ij, is                  ! counter
    real(8), dimension(NS) :: Vnew          ! for updating value function

    ! tolerance level for VFI; maximum number of iterations and speed of updating set in ModuleSolveEqm
    tolV = 1D-10            

    ! Outer loop in a, inner loop in h
    xra_hmax = min(amax,S(:,1)*(1d0+r) + wge*S(:,2)*hmax - TAX(wge*S(:,2)*hmax,r*S(:,1))) ! No unemployment benefit, I choose h
    
    errV = 10d0             ! initialize error value function
    do itV = 1, maxitV

        ! Find optimal a' with golden search
        !$OMP PARALLEL DO DEFAULT(SHARED)
        do is = 1, NS
            a_pol(is) = SolveBell_wa_hopt_sc(amin,xra_hmax(is),is)
        enddo
        !$OMP END PARALLEL DO
        
        ! value associated with optimal choice given by golden
        !$OMP PARALLEL DO DEFAULT(SHARED)
        do is = 1, NS
            Vnew(is)  = VF_eval_wa_hopt_sc(a_pol(is),is,1)  
        enddo
        !$OMP END PARALLEL DO

        ! Check convergence
        errVold = errV                      ! old error value function
        errV    = maxval(abs((V-Vnew)/V))   ! error value function [norm2(V-Vnew)]
        derrV   = abs(errV - errVold)       ! change in error value function

        ! Howard improvement steps
        do ij = 1,120 
            V=Vnew              ! reset value function
            call ComputeVe      ! compute expected value function
            
            ! update value function without optimizing policies
            !$OMP PARALLEL DO DEFAULT(SHARED)
            do is = 1,NS
                Vnew(is) = VF_eval_sc(h_pol(is),a_pol(is),is)
            enddo
            !$OMP END PARALLEL DO
        enddo
        
        if (errV <= tolV) then  ! value function converged           
            !print *, 'Value function converged at itV = ', itV
            exit 
        else                    ! update value function
            V = ww*V + (1.0D0-ww)*Vnew      ! slow updating
            call ComputeVE                  ! compute expected value function
        end if
                
    end do

    end subroutine Solve_GLO_LOC


end module ModuleSolve_GLO_LOC

